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ABSTRACT 




An algorithm to determine the mean pressure field for elliptic flow calculations 
with the probability density function (pdf) method is developed and applied. The 
pdf method is a most promising approach for the computations of turbulent reacting 
flows. Previous computations of elliptic flows with the method have been performed 
in conjunction with conventional finite-volume based calculations that provided the 
mean pressure field. The algorithm developed and described here permits the mean 
pressure field to be determined within the pdf calculations. The pdf method 
incorporating the pressure algorithm is applied to the flow past a backward-facing 
step. The results are in good agreement with data for the reattachment length, 
mean velocities, and turbulence quantities including triple correlations. 


INTRODUCTION 


The understanding and modeling of turbulent reacting flows is of great import- 
ance due to the occurrence of such flows in most practical combustion devices. The 
probability density function (pdf) method offers several advantages over conven- 
tional turbulence models for the computations of turbulent reacting flows (ref. 1). 
In this method a modeled transport equation for the joint pdf of velocities and 
scalars is solved using the Monte Carlo technique. The primary advantages of the 
pdf method are that, in the pdf transport equation, terms representing convection 
and reaction appear in closed form and need not be modeled. Hence, conventional 
turbulence models such as the k-e or second-moment closure models which are based 
on gradient-diffusion assumptions are entirely avoided. Also, any arbitrarily com- 
plex finite-rate reactions that are beyond the reach of conventional turbulent re- 
action models can be treated without approximation or restrictive assumptions. 

Until recently, most of the flows studied using the pdf method have been two- 
dimensional (2-D parabolic or 1-D time-dependent) in nature. These studies (refs. 
2-7 to name a few) have demonstrated the potential advantages of the pdf method es- 
pecially for reacting flows. There are no theoretical limitations for applying the 
pdf method for complex (2-D or 3-D) elliptic flows. However, a suitable algorithm 
for extracting the mean pressure field within the pdf solution algorithm is needed. 

A pdf method for elliptic flows was first demonstrated by Anand, et al . (ref. 

8) for a 2-D flow. The method was applied to 3-D in-cylinder flows by Haworth and 
El Tahry (Refs. 9, 10). In this method, the Monte Carlo (MC) calculations for the 
pdf are combined with conventional finite-volume (FV) calculations (ref. 11) for 
mean fields. The mean pressure field and the turbulent time scale are supplied to 
the MC calculations by the FV calculations. In turn, the MC calculations supply to 
the FV calculations the Reynolds stresses extracted from the pdf solution, thereby 
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avoiding the use of conventional turbulence models. This combined algorithm ex- 
ploits the advantages of both methods to yield an economical algorithm for pdf cal- 
culations of elliptic flows. This method was intended to demonstrate the feasi- 
bility of pdf calculations for such flows. 

However, for more complex flows, especially variable density and reacting flows 
with fast or finite-rate multistep chemistry, the coupling between the two methods 
becomes quite complex and it becomes more difficult to fully exploit the advantages 
of the pdf method. Therefore, it is desirable to perform pdf calculations inde- 
pendently, without recourse to FV calculations. This paper describes an algorithm 
by which pdf calculations are performed independently and economically (in terms of 
computer resources required) for complex elliptic flows. While the present work 
focusses on the determination of the mean pressure field within the pdf calcula- 
tions, parallel work is in progress (ref. 12) that deals with the determination of 
the turbulent time scale. For the present work, the turbulent time scale is sup- 
plied to the calculations as described later in this paper. The relevant equations 
for determining the mean pressure are presented, the solution strategy and the 
numerical scheme are developed and suitably implemented in the pdf calculations. 

The method is applied to a sample elliptic flow, namely the recirculating flow be- 
hind a backward-facing step, and compared against data of reference 13. 

THE PDF METHOD 

A very brief description of the pdf method is presented here. The reader is 
referred to reference 1 and other studies listed in the introduction for more de- 
tails . 

The joint pdf, f(V, X, t), is the probability density of the simultaneous 
event U(x, t) = V and ^(x t) = yi , where U is the velocity vector, | is a set 
of scalars, and V and >Jj are independent variables in the velocity-scalar space. 
Starting from the usual conservation equations for mass (continuity), momentum and 
the scalar quantities, the transport equation for the joint pdf can be derived as 
described in reference 1. In this equation, the terms involving convection, reac- 
tion, body forces and the mean pressure gradient (including the variable density 
effects in those terms) appear in closed form and need not be modeled. Hence, the 
use of conventional turbulence and reaction models is avoided. However, the terms 
in the pdf equation representing the effects of viscous dissipation, fluctuating 
pressure gradient and molecular diffusion need to be modeled. A Lagrangian view- 
point is adopted in modeling and solving the joint pdf equation. The modeled pdf 
transport equation is solved by the Monte Carlo technique. 

In the Monte Carlo solution technique, notional particles representing fluid 
particles are distributed throughout the solution domain. Each particle is attri- 
buted with values for the velocity components, spatial coordinates, and scalar 
quantities. Each of these values evolves according to its corresponding Lagrangian 
evolution equation which incorporates modeled terms where needed. The solution of 
these evolution equations constitutes the solution of the pdf transport equation 
(ref. 1). Means of any functions of the independent variables are determined by 
appropriately summing over the values of all the particles in small spatial sub- 
volumes and fitting splines over the whole domain to these sums. 

In the context of the scope of the present study, the scalar quantities are ig- 
nored in the presentation except to note that the density p is a function of the 
scalar variables. During an interval of time dt, the change in particle position 
x* is given by 
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( 1 ) 


>v >v 

dx. = U.dt, 

1 1 

where U* is the particle velocity. With the use of the simplified Langevin model 
(refs. 2, 3) the increment in U* in the interval dt is given by: 

dlL = - \ ~ ^ dt - (; + 7 C ) (U* - U.>) co dt + (C kco ) 1/2 dW.(t) (2) 

x p* dx. 2 4 o l l o l 

r l 

where angled brackets indicate mean quantities, tilde denotes density-weighted 
^Favre) means, <P>(x*) is the mean pressure, p* is the density of the particle, 

U(x*) is the Eulerian mean velocity, k is the turbulent kinetic energy, co(x*) is 
the turbulent frequency, C Q is a universal constant and W(t) is an isotropic 
Wiener process. The last two terms in Eq 2 jointly represent the effects of viscous 
dissipation and fluctuating pressure gradient. The value C Q = 2.1 was determined 
by Anand and Pope (ref. 2) and has been used in the present study. The quantities 
U and k are easily extracted from the pdf solution at time t. In previous studies, 
the mean pressure gradient has been determined through boundary-layer assumptions 
or has been supplied by the accompanying finite-volume calculations. Its determina- 
tion within the pdf calculations for elliptic flows is the topic of this paper. 

The determination of co is also discussed. 

PRESSURE ALGORITHM 


Overview 


The pressure algorithm is constructed for steady-state, constant- or variable- 
density, high-Reynolds-number flows. The equation for mean pressure is the Poisson 
equation (Eq. 4) derived from the mean momentum equation (Eq. 3): 
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Since the right-hand side of the equation can be evaluated from the pdf solution, 
Eq. 4 can be solved for the mean pressure <P> with appropriate boundary conditions. 
Starting from arbitrary initial conditions, if an iterative or psuedo-time marching 
algorithm is used to reach the steady-state for variable-density flows, then, be- 
cause the transient terms are absent from Eq. 4, the mean pressure given by Eq. 4 
is not correct until the steady-state is achieved. This is true even for constant- 
density cases unless the initial velocity field satisfies the continuity equation 
and Eq. 4 is solved in a coupled manner with the mean momentum equation which is 
not solved directly in the pdf method. A consequence is that the mean continuity 
equation 
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will not be satisfied. Given a field that does not satisfy Eq. 5, a velocity cor- 
rection AU can be obtained by 
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( 6 ) 


AU. 


i <p> 



By requiring that {<pUp> + <p>AUp} be divergence free, we obtain a 
Poisson equation for <f>: 
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(7) 


These observations suggest the following basic algorithm: every so-many steps 


1. solve Eq. 7 for $ 

2. add the velocity correction obtained from Eq. 6 

3. solve Eq. 4 for <P> 

The spirit of the algorithm is similar to that of the SIMPLER algorithm used in 
conventional FV calculations (ref. 11). However, the solution strategy and the im- 
plementation in the pdf calculations are markedly different. 


Boundary Conditions for <P> 

The mean pressure <P> is uniquely determined (to within an arbitrary constant) by 

Eqs . 3 and 4. That is, given Rjj E - <pUpUj>, <P> can be determined from 

the Poisson equation (Eq. 4) with the Neumann boundary conditions provided by Eq. 3 
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where n is the outward pointing normal and n is a coordinate in that direction. 


The numerical solution of a Poisson equation with Neumann conditions everywhere is 
somewhat perilous. Unless suitably treated, the coefficient matrix is singular; 
and a solution exists only if the boundary conditions are precisely consistent with 
the source. To avoid these difficulties, and to make the same boundary conditions 
applicable to $ (see below), a Dirichlet condition is used at the outflow bound- 
ary. There the flow is approximately parallel (to the x-axis) and fully developed, 
and so the lateral momentum equation becomes (approximately) 
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where v" denotes the lateral velocity fluctuation. Hence the Dirichlet condition 
used at the outflow boundary is 


<P> = - <p> v" 2 + P 0 (10) 

The arbitrary constant P Q in the pressure solution is chosen each time so as to 
match a known pressure at any given point in the solution domain. 

Boundary Conditions for <1? 

At walls, planes of symmetry, and inflow boundaries, the normal component of the 
velocity correction should be zero. Hence the Neumann condition 
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is appropriate. At the outflow boundary, specifying the Dirichlet condition, $ = 0, 
results in there being no tangential velocity correction, but allows a correction 
to the normal velocity. This is necessary if the mass flow rate in differs from 
the mass flow rate out. 


NUMERICAL SOLUTION 

In the context that the Monte Carlo calculations are grid-free and the mean 
fields are represented by spline-fits, it is advantageous to solve the Poisson 
equations (Eqs. 4 and 7) in a way that is completely consistent with the spline 
representation of the mean fields. The following discussion is in the context of 
2-D flows, but its extension to 3-D flows is straightforward. 

Let g(x,y) be any function within the solution domain. Then its spline rep- 
resentation is 


M 


N 


_aB 


g(x,y) = l l G^ ,J a (x)b (y), 
a = 1 B = 1 a a 


( 12 ) 


where a a (x)(a = 1, M) are the basis functions in x, bjj(y)(B = 1, N) are 
those in y, and G a ^ are the spline coefficients. Let be the spline 
coefficients of -<pU^>, and let be those of -<pU£U^>. With 

the substitution of the appropriate spline representations, the Poisson equations 
(Eq. 7 and Eq. 4) become 
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and 
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where ( l )a ^ and P a ^ are the spline coefficients of $ and <P>, and primes 
denote differentiation. 


(14) 


Equations 13 and 14, along with the appropriate boundary conditions can each be 
represented by 

A.s = b, (15) 

where s is the solution (spline coefficients of <P> or 3>) and A is the co- 
efficient matrix. Since both <P> and $ satisfy the same types of equations and 
boundary conditions, their coefficient matrices are identical. Further, since A 
is a constant matrix (independent of R-n and Qj[), its LU decomposition need 
be performed once only, and then Eq. (15) can be solved many times by back substi- 
tution. The source vector b is different for <P> and and changes from iteration 
to iteration as the steady-state is approached. 


351 



The structure of matrix 4 is shown in Figure 1, with x's indicating the only 
entries in A when cubic B-splines (refs. 1, 14) are used. A special band solver, 
which exploits the structure of the matrix A was developed to minimize the storage 
and computational requirements. The solver uses LU decomposition and Gaussian 
elimination with partial pivoting. As pointed out earlier, the LU decomposition 
is performed only once and the Poisson equation is solved repeatedly by back sub- 
stitution. The work to perform back substitution is of order MN(4M +5) operations: 
just 34,000 for N = M = 20 used typically (M = N = 21 in the present study). 

The majority of the work is done in calculating b. This involves forming 
splines for <pU^> and <pU^Uj> with appropriate boundary conditions. These 
splines need to be very accurate especially since their derivatives and second 
derivatives are to be evaluated (see Eqs. 4 and 7). Also, the boundary conditions 
used in forming these splines have a significant influence on the solution of the 
Poisson equations. Both the issues need careful consideration. 

In general, the splines in the pdf calculations are formed using the cross-vali- 
dation procedure (refs. 1, 14), by which splines are fit to two independent sets of 
sums (or samples) from each of the spatial bins such that the resulting spline is 
the best fit for both sets of samples. This procedure considerably reduces the er- 
ror in the spline fit with respect to the actual function, compared to the error 
with a simple least-square spline fit with a single set of samples. For spline-fits 
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Figure 1. Structure of the coefficient matrix A (MN x MN) for sample 
values of M and N. There are 5 bands* each with a width of 5 and 
separated by a width of M -5. The total bandwidth is 4M 4- 5. 
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Table 1: Boundary 

conditions for splines needed to 

evaluate b 

Quantity 

q 

Inflow boundary 
normal to x 

Outflow boundary 
normal to x 

Wall 
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subscripts x and y denote differentiation with respect to x and y, respectively 
,oV indicates that no boundary condition is specified 


in two-dimensions, four independent samples are needed. However, even with cross- 
validation, although the first derivatives are acceptable, the second derivatives 
evaluated from the cubic spline fits can have errors significant enough to adversely 
affect the solution of Eq. 4 and hence the convergence of the whole pressure 
algorithm. Hence for Eq. 4, splines are fit for Rij as well as to its appropriate 
first derivative from which the required second derivative is evaluated by a single 
differentiation. From the same raw sums for Rij four different cross-validated 
splines are formed for Rij using four different sets of values for M and N (the 
number of basis functions). The first derivatives evaluated from each of the four 
splines form the samples, although not entirely independent, for the spline-fit for 
the first derivative. This procedure significantly reduces the errors in the second 
derivatives evaluated and results in a stable solution of Eq. 4. The four splines 
were formed in the present study with values of M (=N) equal to 17, 19, 21 and 25. 

It is also necessary to specify a consistent set of boundary conditions for the 
spline fits. The quantities splined for solving Eqs . 4 and 7 and the boundary con- 
ditions used are listed in Table 1. Note that the density has been assumed con- 
stant. The table should be applicable to variable-density flows also, with the 
conventional means replaced by Favre means. Additionally, boundary conditions for 
<p> will be needed — a zero normal-gradient condition at all boundaries is sug- 
gested. 


IMPLEMENTATION 

The development and the implementation of the algorithm to solve for the pres- 
sure field within the Monte Carlo solution algorithm were evolutionary processes. 
Primary concerns were to obtain an accurate pressure solution and also keep the 
computational times low. The algorithm can be summarized as follows: 
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1. form the matrix A and perform its LU decomposition 

2. initially set <P> = 0 everywhere 

3. march the pdf calculations by a step in psuedo-time 

4. perform velocity correction; i.e. 
o form splines of <p> and <pU^> 

o solve Eq. 7 by back substitution for $ 
o correct particle velocities by Eq. 6 

5. on the first step and subsequently, every S steps update pressure; i.e. 
solve for the pressure correction A<P> from Eq. 4 


v 2 a<p> = 9 R ii - v 2 <p> old , 

3x . 3x . 
i 3 


and reset the mean pressure by 


„ new „ old . „ 
<P> = <P> + rA<P>, 


(16) 


(17) 


where r is a relaxation parameter. 

6. if steady-state has not been reached, go to step 3 

7. stop 

Two parameters S and r related to pressure updates, are introduced. For the 
case studied, neither parameter influenced the results significantly. The pressure 
updates are considerably more expensive than velocity corrections. The pressure 
updates need not be performed frequently since the velocity corrections, which cause 
a similar effect as the pressure corrections, are made every step. Also, the 
velocity corrections are quite small from step to step (especially if continuity is 
satisfied initially) so that the pressure updates are not necessary every step. It 
is better to update the pressure after the velocity field has evolved appreciably, 
say every 10 to 20 steps. For the present study, values between 0.5 to 0.8 were 
used for the relaxation factor without any apparent effect on the results or the 
number of steps needed for convergence. 

Although convergence is not difficult to achieve, it is difficult to monitor. 
Several parameters were considered as indicators of convergence. But, due to the 
stochastic nature of the calculations and due to the fact that the solution evolves 
slowly from step to step, a robust and reliable parameter has not been found. In 
the present study; the convergence was inferred by comparing the results at the 
200th step with those at the 220th and the 250th step and concluding that the re- 
sults were unchanged within statistical errors. 


TEST CASE AND INITIAL CONDITIONS 

The pressure algorithm and the pdf calculations were tested for the flow past a 
backward facing step with data from reference 13. The schematic of the flow is 
shown in Figure 2. The step height is H and the centerline mean axial velocity in 
the channel is U re f . The fluid used in the experiments was water, and the 
Reynolds number based on H and U re f was approximately 12,000. 

The inlet velocity pdf was prescribed to be Gaussian for each of the velocity 
components with the means and variances being those of fully developed turbulent 
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H = 0.0762 m W 3 /(W 3 -H)= 1.43 

W 3 = 0.254 m U re f = 0.196 m/s 

TE89-4905 


Figure 2. Schematic of the flow considered. Backward-facing step 
studied experimentally by Pronchick and Kline (1983). 


flow in the inlet channel as calculated by the k-e model. The inlet turbulent 
kinetic energy is distributed into the three components according to k = <u'^> = 
2 <v' 2> = 2<w'^>. The turbulence levels at the inlet have little effect on the 
final results since most of the shear and turbulence production occur in the re- 
circulation zone within the solution domain. The initial pdf in the solution domain 
was set equal to the inlet pdf and the initial mean pressure was zero everywhere. 

The solution domain extended from the step (x/H = 0) to x/H = 15 and between the 
two walls in the y-direction. 

The turbulent frequency (w = e/k, where e is the dissipation rate of k) 
is calculated at each step from the current value of k by 

w = k°> 5 / 2 ., (18) 

where the turbulent length scale 2 . is supplied from the k-e solution of ref. 8 
(see Eq. 19) and is held constant throughout the calculations: 

ft. = k 1 * 5 ^. (19) 

The choice of supplying 2, was preferred over the choices of supplying to or c 
from the k-e solution since the profiles of 2 . were nearly uniform across the 
flow (approximately 0.2 W 3 ) except near the walls where the values decrease. This 
way, the w used in the calculations (Eq. 18) would be more consistent with the 
current values of k. As mentioned earlier, work is in progress (ref. 12) that deals 
with the determination of w within the pdf calculations. 

RESULTS AND DISCUSSION 

The number of particles used in the simulation was 60,000. The number of 
spatial bins used was 1500 (50 in x X 30 in y). The computations were performed 
for 250 steps with pressure updates every 20 steps and the pressure relaxation fac- 
tor of 0.8. The time interval for each step was 0.3 H/U re f. This time interval 
was less than one-tenth of the turbulent time scale (inverse of co) over most of 
the flow domain except for small regions near the wall where it was between one- 
fifth and one-tenth of the turbulent time scale. The calculations reached a steady- 
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state around the 200th step. As noted earlier, a more precise way of determining 
convergence needs to be identified. 

The final results shown here (except for the mean pressure) are calculated from 
spline-fits to sums formed over 50 steps after steady-state has been reached. This 
contributes significantly to reducing the statistical error in the splined results 
especially for the triple and higher correlations which are more prone to such 
error. This procedure, as an alternative to increasing the number of particles in 
the simulation, allows the computer storage requirement to be kept at a modest 
level . 

The computed reattachment length of 6.65H is in good agreement with the experi- 
mentally observed reattachment length (x RE ) of 6.75H. The results for mean 
velocities and various turbulent quantities are shown in Figures 3-11. The re- 
sults are shown for different axial stations normalized with respect to the experi- 
mental reattachment length: 

x* = (x - x R g) /x RE (20) 

Overall, the results are in excellent agreement with data both in magnitude and in 
qualitative trends. It is noteworthy that the turbulence quantities, especially 
the third moments (which cannot be calculated with the k-e model, and are modeled 
in second-moment closures) are in good agreement with data. The slight disagreement 
between some of the computed turbulent quantities and the data in the vicinity of 
y/H=2.0 could be due to the inadequacy of the prescription of the turbulent frequ- 
ency in the present study. It is expected that the use of the stochastic dissipa- 
tion model (ref. 12) will remedy this deficiency. 
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Figure 3. 


Computed mean streamwise velocity profiles 
against data (symbols). 


(lines) compared 
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Figure 4. Computed mean transverse velocity profiles (lines) compared 

against data (symbols). 


x*= -0.71 -0.26 -0.01 0.47 



Figure 5. Computed streamwise velocity variance profiles (lines) compared 

against data (symbols). 




y/H y/H 



Figure 8. Computed <u'^> profiles (lines) compared 
against data (symbols). 


x*= -0.7) -0.26 -0.01 0.47 



103 <v’3>/u^ TE90-966 


Figure 9. Computed <v'3> profiles (lines) compared 
against data (symbols). 
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VALUES 



x/H 


1= -0.141 
2= -0.125 
3= -0.108 
4“ -0.0S2 
5- -0.076 
6= -0.060 
7- -0.044 
8= -0.028 
9- -0.011 
10° 0 . 005 
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Figure 12. Computed contours of mean pressure (normalized by pU^ e £). 


The contours of the calculated mean pressures are shown in Figure 12. The value 
at the right bottom corner of the solution domain is held fixed at zero in the cal- 
culations. The magnitudes and shapes of the contours are comparable to those ob- 
tained from previous calculations with the k-e model and the combined FV/MC pdf 
method (ref. 8). This fact, in conjunction with the accurate prediction of the mean 
velocity field, shows that the pressure algorithm is functioning satisfactorily. 

Approximately 1.3 megawords of storage was required for the computations. The 
total CPU time required for the computations (250 steps) on a CRAY XMP-22 was ap- 
proximately 11 minutes. This reflects an order of magnitude reduction over the es- 
timated time it would have taken with a simple straightforward implementation of 
the algorithm. Nearly 55 percent of the computational time is spent in computations 
related to the pressure algorithm, namely forming sums and splines for pressure up- 
dates and velocity corrections, in spite of limiting the frequency of pressure up- 
dates and improving the efficiency of the spline-fitting routines. However, the 
total computational time of approximately 11 minutes and the storage required are 
still very modest for pdf calculations of an elliptic flow. 

CONCLUSIONS 

An algorithm to determine the mean pressure in steady-state constant or vari- 
able-density elliptic flow calculations with the pdf method has been developed and 
implemented economically. The algorithm and the pdf calculations have been tested 
for a sample recirculating flow; namely the flow behind a backward-facing step. 

The computed results are in excellent agreement with data for the mean velocity 
field and various turbulence quantities including triple corrections. The computed 
reattachment length is in good agreement with the experimental value. These com- 
parisons have served to validate the pressure algorithm and its implementation. 

The pressure algorithm will be validated for variable-density elliptic flows as 
well. Further, with the inclusion of the stochastic dissipation model (ref. 12), 
fully independent pdf calculations of complex elliptic flows will be possible. 
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